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1 Introduction 



Q ; Abstract 

^■^ ■ A model of the heart tissue as a conductive system with two interacting pace- 

makers and a refractory time, is proposed. In the parametric space of the model 
the phase locking areas are investigated in detail. Obtained results allow us to 
predict the behaviour of excitable systems with two pacemakers depending on the 
type and intensity of their interaction and the initial phase. Comparison of the 
' described phenomena with intrinsic pathologies of cardiac rhythms is presented. 
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^ . One of the remarkable examples of excitable media is the cardiac tissue. Because the 

stability of its behaviour is essential for living creatures, investigations of processes 
occurring in the cardiac muscle attract a considerable interest of various scientists. 
Owing to a great complexity and a wide variety of the processes taking place in the 
heart, it can be considered with several points of view. One of them relies upon the 
interpretation of the cardiac tissue as an active conductive system. Then, the cardiac 
rhythms are described on the basis of the dynamical systems theory (see, for example 
[Glass et ai, 2002]). 

The excitation wave in the cardiac tissue originates in the sinoatrial node (SA) 
and spreads rapidly in succession over the right atrium, the left atrium, then the 
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atrioventricular node (AV), bundle of His and Purkinje fibers, and finally to the walls 
of the right and left ventricles. The normal rhythm of the heart is determined by the 
activity of the SA node which is called the leading pacemaker (a source of concentric 
excitation waves) or the first order driver of the rhythm. In addition to the SA node 
cells, the other parts of the cardiac conductive system reveal an automaticity. So, the 
second order driver of the rhythm is located in the AV conjunction. The Purkinje fibers 
are the rhythm driver of the third order. Moreover, for some pathological states of the 
heart arising of ectopic pacemakers is typical. 

Appearance of several excitation sources leads to various disorders in the cardiac 
rhythm, i.e. arrhythmias [Schamorth, 1980; Marriot & Conover, 1983; Zipes & Jalife, 
1985; Winfree, 1987; Glass & Mackey, 1988]. Because arrhythmias are dangerous 
diseases of the heart, their investigations have a great importance. Analysis of the 
complex cardiac rhythms on the basis of their interpretation as chaotic phenomena 
can give a clue to the problem of controllability of the complex cardiac dynamics 
and removing the cardiac tissue to the required regime [Goldberger & Rigney, 1988; 
Goldberger, 1990; Garfinkel et al, 1992]. 

As known, some arrhythmias can be presented by the interaction of spontaneous 
nonlinear sources [Glass et al, 1983; Glass & Mackey, 1988; Kremmydas et al, 1996]. 
In some cases, a model of impulse systems describing certain types of such arrhythmias 
is constructed in the framework of the theory of dynamical systems. This model is rep- 
resented by coupled low-dimensional maps (circle maps) which can exhibit a complex 
dynamics. Such an approach is based on the fact that experiments on periodically 
stimulated cardiac cells are in a close agreement with the dynamics predicted by one- 
dimensional circle maps [Courtemanche et al, 1989]. 

In our investigations, a quite general model of two nonlinear coupled oscillators 
describing certain types of cardiac arrhythmias is constructed. The model turns out 
to be a universal one in the sense that it does not depend on the specific form of inter- 
actions, i.e. on the phase response curve (PRC). The experimentally obtained PRC is 
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approximated by a certain polynomial function with plateau. This plateau describes 
a refractory time when the system does not respond to an external action. Note that 
the refractory time plays an important role for the normal cardiac functioning. For 
example, the refractoriness extends over almost the whole period of the cardiac con- 
traction protecting the myocardium from premature heartbeats caused by an external 
perturbation. The refractoriness provides also the normal sequence of an excitation 
propagation in the heart tissue and the electrical stability of the myocardium [Marriot 
& Conover, 1983; Zipes & Jalife, 1985; Winfree, 1987; Glass & Mackey, 1988]. In the 
proposed model, taking into consideration the refractory time possible areas of the 
phase lockings are investigated. Phenomena of the splitting of resonance tongues and 
superposition of the synchronization areas are found. Using the obtained results we 
can define dynamics of the excitable media with two active pacemakers depending on 
the type and intensity of their interaction and the initial phase difference. Moreover, 
generalizing the principles of our construction one can develop a quite general theory 
of excitable media with interacting pacemakers under external actions. This fact has a 
great practical importance because admits to realize the control of the cardiac rhythms 
by external stimuli. 

2 Heart Tissue as a Dynamical System 

In certain cases cardiac arrhythmias can be described as an interaction of two sponta- 
neously oscillating nonlinear sources. This interaction can be considered as an influence 
of some external periodic perturbation on a nonlinear oscillator (with the constant am- 
plitude and frequency). So, for the description of such a situation it is possible to use 
the well-known circle map [Glass et al, 1983; Bub & Glass, 1994; Kremmydas et al, 
1996]: 

x n+l = x n + f (x n ) (mod 1) , 

where phase difference in oscillators and the function / (x) determines a change 

in the phase after action of stimulus. This function is called a phase response curve 
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(PRC). 

One of the most important characteristics of the circle map is a rotation number p. 
It is defined as follows: 

p = lim . 

n— >oo ji 

For stable phase locking N : M the rotation number is rational, p = N/M. If it is 
irrational the system behaviour is quasiperiodic or chaotic. 

Analysing dynamics of the constructed model based on the circle map, it is necessary 
to find a proper analytical approximation of the experimentally obtained phase response 
curve. This allows us to investigate the basic features of the behaviour of the considered 
system. 

Experiments on the recording of phase shifts have been carried out for a quite large 
number of various systems. First of all, we are interested in the PRC experimentally 
obtained from the research of some cardiac tissues. In [Weidmann, 1961; Jalife & Moe, 
1976] measurements of the cycle durations of the spontaneous beating Purkinje fibres 
after stimulation by short electric current pulses have been performed. The found phase 
response curve is shown in Fig.|I] (dotted lines). Taking into account this experimental 
material, it is possible to make the following general conclusions [Weidmann, 1961; 
Jalife & Moe, 1976; Glass et al, 1986; Glass & Mackey, 1988]: 

• after perturbation the rhythm is usually restored (after some transient time) with 
the same frequency and amplitude, but the phase is shifted; 

• depending on a phase the single input can lead to either increasing or decreasing 
of the period of a perturbed cycle; 

• at some amplitudes of stimulus the obvious breaks appear. 

The basic feature of any approximation of the PRC is the dependence on two 
physical parameters: on the amplitude of stimulus and the input phase. In the ideal 
case the other (so-called "internal") parameters can be reduced to them. 
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Taking into account the polynomial function for the approximation of the PRC, we 
construct a model of two mutually interacting impulse active oscillators. 

3 Analytical Model with a Mutual Influence of Im- 
pulses 

Let us consider the system of two nonlinear interacting oscillators (Fig.|J). Suppose 
that the pulse of the first oscillator with period T\ beats at t n , and the pulse of the 
second oscillator (with period T 2 ) beats at r n . Then the moments of time of the next 
appearance of the impulses are defined as follows: 

t n +l = tn + T\, 
Tn+1 = T n + T 2 . 

Now, taking into account the change in the period of the first oscillator under the 
influence of the second impulse by the value of Ai((r n — t n ) /TiJ, one can get that 
t n+ i = t n + Ti + Ai(^(r n — t n ) /Txj. By the same manner, for the second oscillator 
T n +i = T~ n + T 2 + A 2 ((t n +i — T n ) /T<zj. Dividing these expressions by T\ we arrive at 
the corresponding values for the phases: 

<fn+l = <fn+ T^tAx (5 n - (fi n ) , 
^1 

T-2 1 A / t n T\ 1 . . T n \ 

o n +i = <)„ + — + — A 2 I — + — + —A 1 (d n - tp n ) - — j . 

Here ip n = t n /T\ is a phase of the first perturbed oscillator with respect to the un- 
perturbed one (with period and S n = r n jT\ is the phase of the second perturbed 
oscillator with respect to the same first oscillator with period T\. Using parameters 
a = T 2 /Ti and Ai/Ti = fi, A 2 /T x = f 2 one can write: 

Vn+l = <£n + fl {S n ~ Vn) , 

Sn+1 = $n + a + f 2 (~ ((f n + 1 + fl {S n ~ Vn) ~ S n ) \ . 

Now, omitting intermediate calculations, for the final expression of the phase difference 
in the oscillators we get 

x n+1 = x n + a + f 2 Q (1 + A (a5„) - x n )) - /x (x n ) (mod 1), (1) 



where x n = 5 n - ip n . 

It is obvious, that the PRC changes the form depending on the amplitude of the 
external stimulus. In the simplest case this dependence can be considered as the 
multiplicative relation. Then the phase response curves can be written as 

h = jh (x) , f 2 = eh (x) , 

where h (x) is a periodic function and h (x + 1) = h(x). In such an assumption, the 
expression (|1|) takes the form: 
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n+l 



x n + a + eh (— (1 + jh (x n ) — x n )^ — 7/1 (x n ) (mod 1) 



In the present paper we dwell on the investigations of the map (|2|) with a polynomial 
function h(x). The obtained results are the continuation of our previous works con- 
cerning modeling certain cardiac arrhythmias [Loskutov, 1994; Loskutov et al, 2002]. 

4 Phase Diagrams for Unidirectional Coupling of 
Oscillators 

First of all let us analyze the situation when the permanent inputs act on the nonlinear 
oscillator, i.e. f'2 (x) = or e = 0. As an analytical approximation of the experimental 
curve in Fig.|I|, let us consider the following polynomial function: 

h( x ) =Cx 2 (1-x) 2 . (3) 

The normalizing factor C we choose in such a way that the amplitude of h(x) is equal 
to 1, so that C = 20\/5 (see Fig.|l|, solid line). Then taking into account the refractory 
time 5 the map (Q) we can write as follows: 



x n + a, < x n < 5, (mod 1) 

x n & 
1-5 



Xn+1 ^ Xn + a + C-yhl^-Pl , 5<x N <l, (mod 1), (4) 



where h(-) is determined by (j3|). Now let us compare several cases with different values 
of the refractory time. First of all consider the case without the refractory period, i.e. 
5 = 0. The phase locking regions in the parametric space (a, 7) obtained by numerical 



analysis are shown in Fig.|3|a. Without loss of generality, in this Figure we choose 
a G [1,2]. Different colours define the phase locking areas with the multiplicity N : M, 
where N cycles of external stimulus correspond to M cycles of nonlinear oscillator. 
One can see that "tales" of the main locking regions are slightly splitted and overlap 
each other at large 7. Note that as follows from the analysis of the system (f|) with 
5 = 0.1 (Fig||b), introduction of the refractory time leads to the extension of the phase 
locking areas and significant splitting and overlapping their "tales". 

In Fig.|]a the numerically constructed phase diagram in the case of 5 = 0.3 is 
presented. For the comparison, in the given Figure the same N : M stable phase 
lockings as in Fig. 4 are shown. One can see that when the value of the refractory time 
is growing, the 2 : 3 phase locking area is increasing with simultaneous decreasing of 
the 1 : 1 and 1 : 2 areas. 

The phase locking regions in the case of 5 = 0.5 are shown in Fig.[|b. This phase 
diagram is qualitatively different from pictures considered above. The form of 2 : 3 
phase locking area is stretched and looks like an arrow. The forms of 3 : 4 and 3 : 5 
regions also resemble arrows in the case of 6 = 0.7 (Fig.ga). At 5 = 0.9 all phase 
lockings are degenerated into vertical lines. This situation is presented in Fig.[|b. Note 
that in the case of 5 = 1 (i.e. the system does not respond to the external action) there 
is no any dependence on the stimulus amplitude 7. 

5 Phase Diagrams for Systems with Mutual Inter- 
action 

In this section the system (|2]) at 5 = 0.1 is considered. The analysis is performed in 
the (7, a) and (7, e)-parametric spaces. 

5.1 Phase locking areas in the (7, a)— space 

Let us consider the case of a mutual interaction of two impulse systems. Assume that 
the influence of the first oscillator on the second one is small enough, for example, 
e = 0.1. The corresponding phase diagram displaying the possible behaviour regimes 
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of the system for 5 = 0.1 is shown in Fig.|G]a. One can see that the mutual interaction 
leads to the deformation and the splitting of the phase locking areas. Note that even 
for small values of the amplitude of the second stimulus 7, the overlapping of the main 
phase lockings takes place. Thus, the system dynamics becomes multistable. This 
corresponds to the situation when the limit state of the map depends on an initial 
phase difference x n . The growth of the refractory time in the model with e = 0.1 leads 
to more deep distortion in the forms of the main tongues and disappearance of the 
splitting areas. If, however, we increase the influence of the first oscillator up to, for 
instance e = 0.5, then one can see a very complicated structure with much more deep 
deformation of the main phase locking areas (see Fig.pj). For example, the 1 : 1 area 
is degenerated into a narrow strip, whereas the 1 : 2 phase locking area increases due 
to appearance of long narrow tongues. 

The numerical analysis shows that at the growth of e up to approximately 0.5 the 
occupied by the resonance zones area becomes larger. At the same time, the shape 
of the phase lockings is complexified, and their location is changed. This leads to 
the almost full mixed picture, such that we can find zones of various multiplicity in a 
small neighborhood of almost any point (7, a). However, at the given values of e the 
self-similarly structures are clearly observed. 

Additionally, we have found that at the further growth of the nonlinearity parameter 
e the resonance zones decrease, occupying the less space. In this case the mixing of 
resonance tongues also takes place. Thus, increasing the effect of the influence of 
oscillators leads to the mixing of initially quite regular structure in (7, a)-space. 

5.2 Phase locking regions in the (7, e)— space 

Now we construct the phase diagrams of the interacting oscillators in the space of 
influence amplitudes, i.e. (7, e). In the first instance, let us consider a = 2 (Fig.|7|a). 
This value of the ratio of periods means that for 7 = ^ = the rotation number is 
rational, and the dynamics of the system is periodic with the 1 : 2 phase locking. 
Although at the growth of nonlinearity it is possible to obtain the phase locking areas 



with another multiplicity, even at the large values of 7 and e the system behaviour is 
periodic with the 1 : 2 phase locking. 

Another situation is observed at a = tt/2. Here the rotation number is irrational 
at zero stimulus amplitudes, and the system exhibits the property of quasiperiodicity 
or chaoticity. However, at increasing the nonlinearity the possibility of the appearance 
of the periodic behaviour exists (Fig.[7|b). Here, as for a sufficiently large e, one can 
observe decreasing of the area occupying the resonance zones. Therefore, at irrational 
values of a the probability of the complex behaviour of the system fl2|) is a large enough. 

6 Analogy with Pathological Heart Rhythms 

Summarizing, we shall try to make an analogy between the obtained results and the 
pathological states of the cardiac tissue. Using the developed models it is possible, for 
example, to describe the interaction of the sinus and the ectopic pacemakers, the SA 
and AV nodes and impact of an external perturbation on the sinus pacemaker. 

Similar case of mutual interaction of the SA and AV nodes was investigated in [di 
Bernardo et al, 1998; Signorini et al, 1998]. The authors modeled the AV node as a 
van-der-Pol oscillator and the SA node was considered as a certain modified van-der- 
Pol oscillator. The obtained waveforms satisfactory replicated the action potentials of 
the SA and AV node cells. A bifurcation analysis performed in these works showed a 
possibility to reproduce and classify various types of cardiac pathologies. 

Now let us consider some types of arrhythmias one can predict on the basis of 
our model. If the first pulse oscillator is presented as the SA node and the second 
one is considered as the AV node, then we come to conclusion that some stable phase 
lockings correspond to the cardiac pathologies which are observed in a clinical practice. 
In this case among various obtained lockings one can reveal the normal sinus rhythm 
(1:1 phase locking). In addition, in the diagrams we can see the classical rhythms of 
Wenckebach (N : (N — 1) phase lockings) and N : 1 AV blocks. 

When the first pulse system is considered as the AV node and the second one is 
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presented as the SA node, we obtain the known inverted Wenckebach rhythms found 
in some patients. 

Presence of the wide areas of phase lockings (see Fig.^-0) confirms that in such 
systems it is possible to observe the various kinds of synchronization of two oscillators 
qualitatively corresponding to some types of cardiac arrhythmias. The phase diagram 
allows us to reveal under what conditions of the interaction (i.e. at what values of the 
parameters a, 7, e and 5) one or another type of synchronization exists. Moreover, the 
phase pictures indicate that at the increasing the nonlinearity (i.e. at the growth of 
the parameter 7) areas with various phase lockings are overlapped. The knowledge of 
such regions and the dynamics permits us to remove the system from an undesirable 
mode of synchronization to a more appropriated state by the external action. 

7 Concluding Remarks 

In the presented paper a quite general model of two nonlinear interacting impulse 
oscillatory systems is developed. On the basis of this model it is possible to predict 
some types of cardiac arrhythmias. The constructed model is a universal one in the 
sense that it does not depend on the chosen interaction type, i.e. on the form of 
the phase response curve. Taking into account the refractory time the possible phase 
locking regions of the polynomial maps which describe a nonlinear oscillator under the 
permanent inputs, are investigated. It is found that involving the refractory time leads 
to the extension of the phase locking areas, significant splitting and overlapping their 
"tales". Moreover, the forms of the phase locking stretched and degenerated 

into vertical lines as the refractory time tends to one. 

Detailed analysis of the phase diagram of the system with two mutually interacting 
oscillators in the (7, a)-space shows that besides splitting of the central tongues there 
is an overlapping of the main regions of synchronization which corresponds to various 
types of cardiac arrhythmias. This bistability is observed even for small enough values 
of the interaction. Increasing the refractory time leads to the distortion in the form of 
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the main tongues and disappearance of the splitting areas. For sufficiently large values 
of the interaction we obtain a very complicated picture, where the phase locking areas 
are interwoven with each other. 

In addition, in the constructed model the phase lockings in the space of the stimulus 
amplitudes are observed. It is found that the interacting oscillators can be synchro- 
nized even if the ratio of their periods is irrational (note, that the probability of this 
phenomenon is a quite small). However, in the case without coupling this would cor- 
respond only to the complex dynamics (quasiperiodic or chaotic). 

The obtained results allow us to predict the dynamics of oscillatory systems de- 
pending on the initial phase difference, on the type and the intensity of the interaction. 
Moreover, using the principle of the construction of the model one can develop a quite 
general theory of the interacting oscillators under periodic perturbation. In this case 
the knowledge of the multistability areas can help to stabilize the system dynamics 
and remove the cardiac tissue to the required type of the behaviour. 
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Figure captions 

Fig.l. The phase response curves: the experimental curve (dotted line) and its ana- 
lytical approximation (solid line). The experimentally obtained phase response curve 
shows a dependence of the duration of perturbed cycle (in %) on the phase of input. 
Fig. 2. A construction of the model of two nonlinear interacting oscillators. 
Fig. 3. The phase diagram of the map (|J): a) 5 = 0; b) 5 = 0.1. 
Fig. 4. The phase locking areas of the map (|4]): a) 5 = 0.3; b) 5 = 0.5. 
Fig. 5. The phase diagram of the map (|j): a) 5 = 0.7; b) 5 = 0.9. 
Fig. 6. The phase locking regions of the system of two mutual interacting oscillators 
with 5 = 0.1: a) e = 0.1; b) e = 0.5. 

Fig. 7. The phase lockings in the space of stimulus amplitudes (5 = 0.1): a) a = 7r/2; 
b) a = 2. 
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Figure 2: A. Loskutov, S. Rybalko & E. Zhuchkova 




Figure 3: A. Loskutov, S. Rybalko & E. Zhuchkova 
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Figure 5: A. Loskutov, S. Rybalko & E. Zhuchkova 
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